function [groups, C, frames, descriptors, nbx] = test2(imgpath, saveimg)
% Input argument:
%   IMGPATH : Path to image file.
%   SAVEIMG : 1 to save output image in each step.
% Output argument:
%   GROUPS
%   C
%   FRAMES
%   DESCRIPTORS
%   NBX
close all;
if nargin < 1
    imgpath = '/image/S/jpg/Tp_S_NNN_S_B_art00015_art00015_20062.jpg';
end
if nargin < 2
    saveimg = 0;
end
img = imread(imgpath);
gray = rgb2gray(img);
scaled = double(gray / 256);
%% SIFT Features
hold off;
imshow(img);
disp('Step 1: SIFT');
[frames,descriptors] = sift(scaled);
hold on;
plotsiftdescriptor(descriptors, frames);
if saveimg == 1
    print -djpeg step1
end
pause(0.5);
disp('Step 1: SIFT finished');

%% Match Features
disp('Step 2: GETNEARESTNEIGHBOURRATIO');
nb = getNearestNeighbourRatio(frames, descriptors);
nbx = nb;
for i = 1:length(nb)
    if nbx(i) <= i
        nbx(i) = NaN;
    end
end
nbs = [1:length(nbx);nbx];
nbs = nbs(:, ~isnan(nbs(2,:)));
hold off;
plotmatches2(img, frames, nbs);
fprintf(1,'Total Matched keypoints: %d\n',size(nbs,2));
if saveimg == 1
    print -djpeg step2
end
pause(0.5);
disp('Step 2: GETNEARESTNEIGHBOURRATIO finished');

%% Cluster
disp('Step 3: CLUSTER');
Z = linkage(frames', 'average', 'euclidean');
%       'single'    --- nearest distance (default)
%       'complete'  --- furthest distance
%       'average'   --- unweighted average distance (UPGMA) (also known as group average)
%       'weighted'  --- weighted average distance (WPGMA)
%       'centroid'  --- unweighted center of mass distance (UPGMC)
%       'median'    --- weighted center of mass distance (WPGMC)
%       'ward'      --- inner squared distance (min variance algorithm)
maxClust = 10;
if maxClust > size(nbs, 2) / 3
    maxClust = ceil(size(nbs, 2) / 3);
end
C = cluster(Z, 'maxclust', maxClust);
hold off;
imshow(img);
hold on;
nCluster = max(C)
countCluster = histc(C,[1:nCluster]);
availCluster = countCluster > 3; % More then 3 points to fit a line
for i = 1:nCluster
    color = rand(1, 3);
    if ~availCluster(C(i))
        continue
    end
    q = find(C == i);
    p = nbx(q);
    ind = find(~isnan(p));
    p = p(ind);
    q = q(ind)';

    nPoint = length(q);
    for j = 1:nPoint
        plot([frames(1,q(j)), frames(1,p(j))],[frames(2,q(j)), frames(2,p(j))],'--s','MarkerSize', 4, 'MarkerFaceColor', color, 'Color', color);
    end
end
if saveimg == 1
    print -djpeg step3
end
pause(0.5);
disp('Step 3: CLUSTER finished');

%% RANSAC
disp('Step 4: RANSAC finished');
nGroup = 0;
groups = [];

for i = 1:nCluster
    q = find(C == i);
    p = nbx(q);
    ind = find(~isnan(p));
    p = p(ind);
    q = q(ind)';

    if length(p) < 3
        %Minimal sample set dimension = 3
        continue;
    end
    X1 = frames(1:2, q);
    X2 = frames(1:2, p);
    
    [sx, sy, theta, H, CS] = RANSACaffine(X1, X2);
    if isnan(theta)
        % no result
        continue;
    end
    X1 = frames(:, q);
    X2 = frames(:, p);
    X1 = X1(:,CS ~= 0);
    X2 = X2(:,CS ~= 0);
    gb = norm([H(1,3),H(2,3)]);
    
    if (abs(theta) < 0.05) && (gb > norm([32 32])) % no rotate and some shift
        foundSimilar = 0;
        for j = 1:nGroup
            gx = (groups{j}.sx - sx) / (groups{j}.sx + eps);
            gy = (groups{j}.sy - sy) / (groups{j}.sy + eps);
            gt = norm([groups{j}.H(1,3) - H(1,3), groups{j}.H(2,3) - H(2,3)]);
            if (abs(gx) < 0.1) && (abs(gy) < 0.1) && (gt < 10)
                    % same group?
                groups{j}.C = [groups{j}.C, i];
                groups{nGroup}.FrameP = [groups{nGroup}.FrameP, X1];
                groups{nGroup}.FrameQ = [groups{nGroup}.FrameQ, X2];
                foundSimilar = 1;
                break;
            end
        end
        if foundSimilar == 0
            nGroup = nGroup + 1;
            groups{nGroup}.C = [i];
            groups{nGroup}.sx = sx;
            groups{nGroup}.sy = sy;
            groups{nGroup}.H = H;
            groups{nGroup}.theta = theta;
            groups{nGroup}.FrameP = X1;
            groups{nGroup}.FrameQ = X2;
        end
    end
end
fprintf(1, 'Total Group Number: %d\n',nGroup);
disp('Step 4: RANSAC finished');
pause;
%% draw groups
if nGroup < 1 %NO GROUPING RESULT, DO NOT REDRAW LAST PICUTRE
    return
end
hold off;
imshow(img);
hold on;
for i = 1:nGroup
    color = rand(1,3);
    framesp = groups{i}.FrameP;
    framesq = groups{i}.FrameQ;
    px = min(framesp(1,:) - framesp(3,:));
    py = min(framesp(2,:) - framesp(3,:));
    ph = max(framesp(1,:) + framesp(3,:)) - px;
    pw = max(framesp(2,:) + framesp(3,:)) - py;
    rectangle('Position',[px,py,ph,pw],'LineWidth',2,'LineStyle','--','EdgeColor',color);
    qx = min(framesq(1,:) - framesq(3,:));
    qy = min(framesq(2,:) - framesq(3,:));
    qh = max(framesq(1,:) - framesq(3,:)) - qx;
    qw = max(framesq(2,:) - framesq(3,:)) - qy;
    rectangle('Position',[qx,qy,qh,qw],'LineWidth',2,'LineStyle','--','EdgeColor',color);
    fprintf(1,'Matched Group %d: Scale X=%0.3f, Y=%0.3f\n',i,groups{i}.sx, groups{i}.sy);
end
if saveimg == 1
    print -djpeg step4
end
pause(0.5);
end

